Skip to content

Genome-wide TWAS fixes: conversion edge cases, noResult for empty LD blocks, optional gene unit and CV prefilter - #1439

Open
jaempawi wants to merge 5 commits into
StatFunGen:mainfrom
jaempawi:twas-genomewide-fixes
Open

jaempawi wants to merge 5 commits into
StatFunGen:mainfrom
jaempawi:twas-genomewide-fixes

Conversation

@jaempawi

@jaempawi jaempawi commented Nov 10, 2026 •

Copy link
Copy Markdown
Contributor

Five commits needed to run genome-wide TWAS with the converted FunGen-xQTL weights on pecotmr 0.9.6. They have been used for the MSBB eQTL (17,297 genes) and KNIGHT pQTL runs on the 2026 AD GWAS.

  1. twas_weights_conversion.R: handle legacy weights stored for the QC'd variant subset only, and genes with no usable weights, instead of failing the conversion.
  2. twas.R, twas_manifest.R: gwas_sumstats_construct.R writes no GwasSumStats file for an LD block without GWAS variants, so twas.R failed on readRDS and the whole region errored. Genes there now get an empty result with a noResult attribute.
  3. twas_manifest.R, twas_ctwas.ipynb: new --unit gene (one unit per gene spanning its weight window, so GWAS and LD are read per gene) next to the existing region grouping. New optional CV prefilter (--min-cv-rsq, --max-cv-pval) that skips genes where no context/method passes, so their GWAS input is never built.
  4. Hardening and tests:
    • Prefilter only runs when --min-cv-rsq > 0. The CV metric names are configurable (--rsq-option, --rsq-pval-option). Methods without CV metrics (e.g. ensemble) are skipped. Genes with no CV results are kept; genes whose weight file cannot be read are kept with a message.
    • --unit gene stops if a gene has more than one xQTL meta row, and makes the unit ID file-name safe.
    • Conversion checks for all-zero/NA weights before variant ID resolution.
    • New tests: gene-unit manifest rows, repeated-gene error, prefilter keeps an unreadable file, prefilter drops failing genes, and a gene-unit TWAS run on the protocol example.
  5. Edge cases and docs:
    • --unit gene stops early when a gene has no start/end in the xQTL meta, and when two gene IDs would give the same file name.
    • twas.R: a missing GwasSumStats file still means an empty window, but a missing folder now stops as a wrong path. The twas step removes an old GwasSumStats file before rebuilding it, and makes the gene part of output names file-name safe.
    • The notebook Input and Command Interface sections list --twas_unit, --twas_prefilter, --rsq_option and --rsq_pval_option.
    • New tests: prefilter p-value cutoff (keep and drop), missing gene window, file-name safe IDs and clashes, missing GwasSumStats file and folder, and conversion of weights without names, with unmatched IDs and with all-zero weights.

No default changes for existing users: twas_ctwas.ipynb keeps the region unit and the prefilter off. Opt in with --twas_unit gene (outputs are then named <name>.<gene>.<gene>.twas.rds) and --twas_prefilter. Use the prefilter for TWAS-only runs, since the cTWAS steps read the same manifest and cTWAS should use every gene.

Tested on SCC with pecotmr 0.9.6: MSBB eQTL TWAS on all 22 chromosomes gives the same per-chromosome gene counts as the 0.9.5 run, and the TWAS/cTWAS integration and script tests pass (27 tests).

gwas_sumstats_construct.R writes no GwasSumStats file for a block without GWAS variants, so twas.R failed on readRDS and the whole region errored. Such genes now get an empty result with a noResult attribute.
twas_manifest.R gains --unit gene (one unit per gene spanning its weight window, so GWAS and LD are read per gene) next to the existing region grouping, which stays the script default. An optional CV prefilter (--min-cv-rsq / --max-cv-pval) skips genes where no context/method passes, so their GWAS input is never built; weights without per-row CV results are kept for pecotmr to decide. twas_ctwas.ipynb uses the gene unit and the prefilter by default (--twas_unit, --twas_prefilter).
- twas_ctwas.ipynb keeps region as the default TWAS unit, so existing commands give the same outputs; the gene unit is opt-in with --twas_unit gene.
- CV prefilter: off by default in the notebook (--twas_prefilter turns it on for twas runs), since the ctwas steps read the same manifest and cTWAS uses every gene with weights. It now reads the per-method CV results of the pecotmr 0.9.5+ TwasWeights layout (it previously found no metrics there and kept every gene), ignores methods without a CV r2 such as ensemble, uses the twas step's --rsq-option and --rsq-pval-option, keeps genes with no CV metrics, no cvResult column or an unreadable weight file (reported on stderr), treats --min-cv-rsq <= 0 as off, and loads pecotmr before reading the weights.
- Gene unit: make region_id file-name safe and stop when the xQTL meta repeats a gene, since both would collide on output paths.
- twas_weights_conversion.R: skip all-zero/NA weights before resolving variant IDs, and prefer preset_variants_result IDs when their length fits.
- Tests: the notebook twas test pins the region unit; new tests cover the gene unit (same result as the region unit on the toy gene), gene-unit manifest rows, repeated genes, and the prefilter keep/drop rules.
- twas_manifest.R: --unit gene stops early when a gene has no start/end in the xQTL meta, checks for repeated genes before making IDs file-name safe, and stops when two genes would share a file name. The CV prefilter skips cvResult entries that are not lists instead of reporting the file as unreadable.
- twas.R: a missing GwasSumStats file still means an empty window, but a missing folder now stops as a wrong path; the noResult text no longer assumes an LD block.
- twas_ctwas.ipynb: the twas step removes an old GwasSumStats file before rebuilding it, makes the gene part of output names file-name safe, and the Input and Command Interface sections list --twas_unit, --twas_prefilter, --rsq_option and --rsq_pval_option.
- twas_weights_conversion.R: report genes without usable weights with message(), since the warning was lost on quit().
- Tests: prefilter p-value cutoff (keep and drop, with the pval fall-back and ensemble skip), NA gene window, file-name safe IDs and clashes, missing GwasSumStats file and folder, and conversion of weights without names, with unmatched IDs and with all-zero weights.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant